////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  ----------------------------- Ebrahim Foulaadvand, 23 June 2012 ------------------------------------- //
//                                                                                                        //
//  The programme "Oscillator" evaluates the temporal evolution of a simple harmonic oscillator position. //
//  It also computed the oscillator energy at eacg time step and compares it to the anlytic solutiom.     //
//                                                                                                        //
//                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;


double a (double &x); // a is the accelration function

double  tau=0.01,T=20,ac,k=9.,m=1.,x0=1.,v0=0.,EEU,EEC,E0,delEEU,delEEC;

// tau=time step, k=spring constant, ac=instantaneous acceleration, m=oscillator mass, x0=initial position,
// v0=initial velocity, T=simulation time, E=total energy, delE=difference between analytic energy and
// computed one at each time step. EU stand for Euler algorithm and EC stands for Euler-Cromer algorithm.

int N=(T/tau); // N=number of time steps

main(){

vector<double> vEU(N+1,0),xEU(N+1,0),vEC(N+1,0),xEC(N+1,0);

// Arrays x and v store the oscillator's position and velocity time step t. EU stand for Euler algorithm
// and EC stands for Euler-Cromer algorithm.

E0=0.5*k*x0*x0 + 0.5*m*v0*v0; // initial energy assinement.

// ----------- initial conditions --------------------

xEU[0]=x0; vEU[0]=v0; xEC[0]=x0; vEC[0]=v0;

ofstream filexEU("xEU tau=0.01.plt"); // output file for position-time curve in Euler algorithm.
ofstream filevEU("vEU.plt"); // output file for velocity-time curve in Euler algorithm.
ofstream filedeltaEU("deltaEU tau=0.01.plt");// output file for energy difference-time curve in Euler algorithm.

ofstream filexEC("xEC tau=0.01.plt"); // output file for position-time curve in Euler-Cromer algorithm.
ofstream filevEC("vEC.plt"); // output file for velocity-time curve in Euler-Cromer algorithm.
ofstream filedeltaEC("deltaEC tau=0.01.plt"); // output file for energy difference-time vs time in Euler-Cromer algorithm.


//------------------------ Euler Method------------------------

 for(int t=0;t<N;t++){

 vEU[t+1]=vEU[t] + tau*a(xEU[t]);

 xEU[t+1]=xEU[t] + tau*vEU[t];

 EEU=0.5*m*vEU[t]*vEU[t] + 0.5*k*xEU[t]*xEU[t];

 delEEU=EEU-E0;

 filexEU<<t*tau<<" "<<xEU[t]<<endl;
 filevEU<<t*tau<<" "<<vEU[t]<<endl;
 filedeltaEU<<t*tau<<"  "<<delEEU<<endl;

 }//end t

//------------------------ Euler-Cromer Method------------------------

 for(int t=0;t<N;t++){

 vEC[t+1]=vEC[t] + tau*a(xEC[t]);

 xEC[t+1]=xEC[t] + tau*vEC[t+1];

 EEC=0.5*m*vEC[t]*vEC[t] + 0.5*k*xEC[t]*xEC[t];

 delEEC=EEC-E0;

 filexEC<<t*tau<<" "<<xEC[t]<<endl;
 filevEC<<t*tau<<" "<<vEC[t]<<endl;
 filedeltaEC<<t*tau<<"  "<<delEEC<<endl;

 }//end t


 cout<<"end."<<endl;getch();

 return 0;

       }

//----------------------------------------------------------------

 double a (double &x) { return ac=-(k*x)/m; }